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I. MOTIVATIONS 

Although spin glasses have been studied for over 
two decades, many fundamental aspects of these sys- 
tems [1113 remain unclear. In fact, some of the most basic 
questions relating to their equilibrium properties remain 
open. Roughly, there are two schools of thought about 
spin-glass theories. In the droplet /scaling picture 0, , 
at low temperature T there are just two equilibrium pure 
phases, related by the up-down symmetry; there is no 
replica symmetry breaking (RSB) and connected corre- 
lation functions decay to zero at long distance. On the 
contrary, in the mean- field picture 0^|HIEI3)I3j at low 
T there are many pure states that are not related by the 
up-down symmetry; because of this, generic connected 
correlation functions do not decay to zero at large dis- 
tances. 

Much of the theoretical and numerical work on the low 
temperature phase of spin glasses has been focused on the 
presence or absence of replica symmetry breaking. The 
study of correlation functions has not been pushed as 
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much: at best, one of the spin-spin correlation functions 
has been measured :9j and was found to decrease with a 
power law. Nevertheless, there are a number of different 
theoretical predictions. Within the scaling and droplet 
pictures 0, 0, , 2-spin connected correlation functions 
generically decay to zero as , where r is the distance 
between the spins and i/ = —1/6 is the usual thermal 
exponent. For large dimensions d one has 6 « and 
numerical studies indicate that d{d = 3) ~ 0.20. 

The predictions of the mean-field picture have been 
more difficult to obtain: it is necessary to go beyond the 
Sherrington-Kirkpatrick ,llj model because in that sys- 
tem one cannot define distances. Generally one does this 
using replica field theory which allows computations [l^ 
in dimensions d > 6. Note that this approach automati- 
cally forces one to consider disorder-averaged quantities. 
Within this formalism, one finds that 2-spin connected 
correlation functions do not go to zero at large distances 
(because of RSB) and that the large r behavior is of the 
form A + B r^^'^. 

In this work we reconsider these issues theoretically 
and numerically in d = 3. From replica field theory, we 
find the remarkable property that a certain linear com- 
bination of correlation functions decays to zero at large 
distance; furthermore, this decay is "fast" in an interme- 
diate region, as corroborated by our numerical measure- 
ments. 

The outline of this paper is as follows. In Section ^ 
we specify our model and we define the correlation func- 
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tions that we investigate. In Section UTTl we discuss the 
replica field theory approach and we compute the corre- 
lation functions of interest in the tree approximation (no 
loops). We then move on to the numerical simulations. 
The corresponding methods are presented in Section llVI 
while results are presented in Section We end with 
our conclusions. 



II. GENERAL FRAMEWORK 

A. The Model 

At a microscopic level, we consider a c?-dimensional lat- 
tice of Ising spins, Si = ±1, coupled by nearest-neighbor 
interactions. The corresponding Hamiltonian is 

Hj = — ^ JijSiSj — B ^ Si , (1) 

where we have introduced an external magnetic field B. 
Ferromagnetic couplings correspond to > 0, anti- 
ferromagnetic ones to Jij < 0. These couplings are 
quenched independent random variables symmetrically 
distributed around 0. To simplify the analytic part of our 
study, it is most appropriate to take them to be Gaus- 
sian of zero mean, while for our c? = 3 numerical study, 
we use Jij = ±1 for simulational reasons (the difference 
is irrelevant in the temperature region we explore) . 

All the issues we consider in the following concern tem- 
peratures T low enough so that one is within the frozen 
spin-glass phase. Furthermore, we shall always take the 
B limit. Within the scaling/droplet picture, this 
limit leaves one with a single equilibrium phase; on the 
contrary, in the mean-field picture, _B — > simply selects 
one partner for each pair of _B = equilibrium phases 
(because of RSB, there is an infinite number of these 
pairs). 



-B = 0, by symmetry we have (5*^) = 0; that is why we 
have introduced the infinitesimal magnetic field. Note 
that (•) does not correspond to an average inside a pure 
state; (Si)"^ is thus not the Edwards- Anderson order 
parameter. 

In our study, we focus on disorder averages, averag- 
ing our observables over samples where the couplings Jij 
are randomly chosen with their a-priori distribution. We 
thus define 



Gi{i,j) = Ti{i,j) G2(*,j) = r2(i,j) , (4) 

where the overline denotes this disorder average. 

What do we expect for these quantities? Clearly, in 
the droplet/scahng picture, the hmit B ^ leads to a 
single pure phase and so connected correlation functions 
go to zero. Furthermore, the large distance behavior is 
controlled by a zero temperature disordered fixed point. 
Roughly, this leads to effective couplings of "block" spins 
at large scales that have a broad distribution: typical val- 
ues scale as for blocks at distance r. Since this distri- 
bution is bell shaped with a non-zero value at the origin, 
droplets of characteristic scale r are thermally excited 
with a probability going as kTr^^. This feature leads to 
connected correlation functions decaying as . 

In the mean-field picture the low T phase undergoes 
RSB: connected correlation functions do not go to zero 
at large r. Interestingly, the r — > oo limit of the disorder 
averaged (SiSj^ can be expressed simply in terms of a 
moment of configurational overlaps. One defines the spin 
overlap q of two configurations (1) and (2) as 

o(l)o(2) 

^(1,2) ^ ^.^^^ , (5) 

where N is the total number of spins. A straightforward 
computation of {q'^) gives 



B. 2-Spin Correlation Functions 

For any given sample, consider the standard spin-spin 
connected correlation function: 

(S^S,) ~ {S,){S,) . 

This quantity is not gauge invariant: it depends on the 
(arbitrary) choices of z axis orientations at sites i and j. 
To restore gauge invariance, one should consider either 
the square of this correlation function 

T,it,j) = [{S,S,)-{S,){S,)r , (2) 
or the difference of the squares, 

^,i^,J) = {s,s,r-{s,ns,r ■ (3) 

In all these quantities, (•) refers to the (thermal) ensemble 
average with respect to the usual Boltzmann measure. If 



(1) = (6) 

The large volume limit of this quantity can be obtained 
by replacing the correlation functions on the right hand 
side by their long distance limit: in this way one sees that 
the long distance limit of the first contribution to G2 is 
{q^) ■ As we shall see later, the mean-field formalism also 
predicts a very non-trivial relation between the limiting 
values of Gi , G2 and a certain combination of moments 
of overlaps. 

A more subtle question for the mean-field picture is 
how the correlation functions approach their large r lim- 
its. The framework for addressing this question is replica 
field theory, which leads to the conclusion that the 
decay should follow a power law in r. Our goal here is 
to reconsider that computation, focusing in particular on 
linear combinations of Gi and G2. 
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III. REPLICA FIELD THEORY FOR 
GENERALIZED CORRELATION FUNCTIONS 

A. Starting Point 

We consider the linear combination of correlation func- 
tions aGi + (1 — a)G2 which can be written as 



X"(z,j) = {S,S,)^-2a{S,S,){S,){S,) 
+ (2a-l)(^mF, 



(7) 



where a is a parameter that can take real values. The 
brackets indicate thermodynamical averages in the pres- 
ence of the infinitesimal magnetic field, while the overline 
indicates an average over the disorder distribution. In the 
case a = 1, after summing over j, we obtain the usual 
expression for the spin-glass susceptibility x- 

The starting point of replica field theory is the (trun- 
cated) effective Lagrangian describing the replicated the- 
ory Ea 



C = 



a,b p 



6VN 



^ ^ $'^''(Pi)<i>''^(P2)<i>^''(P3) 



(8) 



a,6,c {pi} 



+ E E 1''^''(Pl)*°'(P2)'i>"'(P3)$"''(P4) 



a, 6 {pi} 



where all the sums over momenta are constrained in such 
a way that the total momentum adds to zero. Note that 
to reduce the complexity of the presentation, we have 
used discrete sums over momenta rather than the inte- 
grals that arise in the thermodynamic limit, N oo. 
Expanding around the mean-field solution 

^''\p) = Sp,oqab + r\p) , (9) 

one finds 

(10) 



a<b,c<d p 



where 

^mf(q) ^ N{- 2_^q^t, + - 2^ qabqbcqca + Y2 



ab 



a^b,c 



a.b 



(11) 

M is the matrix of the Gaussian fluctuations and Cint 
is an interaction term including contributions up to the 
fourth order in (j) (see Eq. ©). 

In the context of this replicated field theory, the spin 
correlation function that appears in Eq. |(7J) can be ex- 
pressed in terms of field correlations. For example, we 
have 



(12) 



where {■■.)c is an average performed using the effective 
Lagrangian of Eq. © and where the superscripts on the 
spins are replica indices. In the low temperature phase, 
where one assumes that many states may exist, the ex- 
pression for x"(*; j) leads one to consider all the possible 
ways in which 4 replica indices can appear. Then we 
have the following general expression (in position space, 
or, equivalently in momentum space): 



n{n 



1 ^|^a6;a. 



b 2a ^ ^^f, 



(n-2)(n-3) 



(n-2) 



a,b c^b 

(2a -1) 



(13) 



c.d^a^b 



where we have dropped the dependence of x" and 
of the Cs to lighten the notation. 

We can now separate, as in Eq. jnj, the fluctuating 
part from the mean- field part: 



= qabqcd + C^'-'^^l-j) 
+ loop corrections , 



t<t>f) 



(14) 



where G = M^^ is the free propagator (see Eq. l^lQ\i ). 
and where we schematically denote hy i — j the distance 
between i and j. From this we get 



Xmf + Xjl — 



2 - a 



3n(n — 1) 

ab 



^ qlb 



2^ qabqac\ 
abc 



1 



{G 



ab-.ab 



2a 



nin — 1) -^-^ 

(2a - 1) 
(n-2)(n-3) 

^ ' c,d=jta,b 

loop corrections . 



^ ' c^b 



\ ^ Qab\ac 



(15) 



Note that in the mean-field contribution Xmf sums 
run over all replica indices but we have included only the 
terms which have a finite limit when n ^ 0. 

We now compute the generalized correlation function 
X" considering the mean-field solution together with the 
Gaussian fluctuations around it, disregarding the loop 
corrections in Eq. (|15|) . 



B. Mean-Field Solution and Free Propagators 

The mean-field solution qab is obtained by minimizing 
the Lagrangian Cmf it is thus given by the saddle point 
equation 



2u 

2Tqab + w[q^)ab + -^qlb 







(16) 



This equation is solved using an ultrametric Ansatz for 
the form of the matrix q [a S 111- According to this 
Ansatz the matrix q is hierarchically divided into blocks 
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by means of a procedure based on R steps: at each step 
r the size of the blocks is set to Pr, with Po = n and 
Pr+1 = 1. If two rephca indices belong to the same block 
of size Pr but to two distinct blocks of size Pr+i then qat = 
Qr , and we say that their "co-distance" a H 6 in replica 
space is r. In this way the matrix q is parameterized 
by the i?-step function qr, where r G [0, i?]; note that 
each qr has a multiplicity 5r = Pr — Pr+i- One can also 
represent this procedure with blocks via a hierarchical 
tree, having its root at r = 0, its leaves at r = i?, and 
having p^/pr+i branches at nodes of height r. Replica 
indices lie on the leaves and their overlap is qr if they 
join at a node of height r on the tree. At the end one 
takes the limit R oo and n 0: now the values pr 
tend to a monotonic increasing function that takes values 
in [0, 1]. In this limit, setting for example pr = x with 
X € [0,1], qr becomes a continuous function q{x) and 
we indicate the replica co-distance as a (1 b — x. The 
well-known mean-field solution of Eq. then reads: 



q{x) 
qix) 



w 

2m"' 



X < Xi 



qi = —xi xi < X <1 
2u 




with 



wqi 



ml 



. 



(17) 



FIG. 1: Topologies associated with the replicon sector (top) 
and the longitudinal anomalous sector (bottom). 

divergences (see 0,^3 ^"^^ ^^e Appendix El of this pa- 
per), enabling one to block-diagonalize the fluctuation 
matrix. In this way the propagators G are expressed in 
terms of Replica Fourier transformed variables, whose ex- 
plicit form in terms of the mean-field solution qab can be 
computed. In the next section we will use these expres- 
sions to compute the generalized correlation functions. 



To calculate the free propagator G, one has to evaluate 
the fluctuation matrix 



The Generalized Correlation Functions 



M 



ab:cd 



dqcd 



lab 



(18) 



given the mean-field solution Eq. H17|) . and then invert 
it. To do this, it is convenient to parametrize both the 
fluctuation matrix M and the propagator G = in a 
way that exploits the ultrametric structure of the mean- 
field solution. Each element of M and G is in principle 
labeled by four replica indices. However, due to ultra- 
metricity, three replica co-distances are enough to char- 
acterize the geometrical structure of the four indices. In 
terms of these co-distances we can then identify two pos- 
sible "sectors" for the propagator G"''^'^'^ (and for M"''''^'*): 

i) The replicon sector: when ar\h — cf] d = r 



with u = max(a H c, a H d) , v = max(6 f] c,b H d) 
and u,v > r. 

ii) The longitudinal-anomalous sector: when aClb = r 
and cC] d = s 



G 



ab:cd 



^g: 



(20) 



with t — max(a (1 c,a d,b c,b (1 d). 



These two sectors are illustrated in Fig. 

The inversion of the fluctuation matrix is not simple 
[l^ . It can be carried out using the Replica Fourier trans- 
form when d > 6, since in this case there are no infra-red 



The expression for the mean-field part Xmf easily 
obtained using Eq. H17|l: 



Xmf 



2 - a 
3 

2 - a k;2 
3 4m2 



dx q^{x) — I / dx q{x) 
\Jo 



(21) 



Note that this expression gives for a = 2 a mean-field 
value equal to zero. As we will discuss later, this fact is 
the consequence of a "sum rule" . However, as we shall 
see in the next sub-section, the case a = 2 is peculiar 
also at the level of fluctuations. 

In the Parisi solution, x{q), the inverse of q{x)^ has 
a simple expression in terms of P{q): 



dx(q) 
dq 



where 



This makes clear that 
2-a 



Xmf 



— a 



Note that in the present computation the overlaps are 
constrained to positive values because of the inflnitesi- 
mal positive magnetic field; however, that will not be 
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case for the numerical data which were collected in zero 
field (see the next section). For this reason and to use 
a common notation throughout the paper, we indicated 
the variance of q with respect to P{q) in the presence of 
an infinitesimal field as a\q\. 

For the fluctuating part, with the parameterization 
given in the previous section, and exploiting the defini- 
tion of Replica Fourier transform, we obtain in the i?-step 
Ansatz: 



loo 

a ) 

lySr G!j^ +l.ySrySs G^A 

3 1 ^ r+l,r+l 2 ^ ^ Of 

I J 



where the indicates that we have used an operation 
of Replica Fourier transform with respect to the hatted 
variable, and rG^^^ stands for a pure replicon propagator 
(where the longitudinal-anomalous contribution has been 
taken away, see Appendix 



1. The Special Case a = 2 

In the case of a = 2 the first two lines of terms 
in Eq. H22I) disappear and the calculation simplifies a 
lot. We start with the general expression for the replica 
Fourier transform of the replicon propagator (see Ap- 
pendix : 



(23) 



and thus the first sum gives 



R R ^ 



(Po - 1)4 A - (24) 



The second sum is less straightforward. When we go 
to the continuous limit, the variables Pr,Ps with r, s < 
R assume values in the continuum interval [0, Xi], while 
Pr^i = 1; because of that, it is convenient to separate 
the sum into two different pieces before taking the limit 



of continuous replica symmetry breaking (i? — > cx)): 
R R 



Sr Ss Gy 



\ysry5sG^^= y 

r<s<R-l 

{PR - PR+i) y Sr Gg + {pR - PR+if Gg 



R.R 




r<R-l 

dx / dy Gl^y 



(1 











xi 



dx G^ 



{l-x,f G- 



(25) 



Now we can use the expression of the longitudinal Fourier 
transform obtained in Il2l: 



^6 



sinh ■ 



cosh 



where 
Eq. 123) 



p^p cosh ^ 
x\-V 1 - xi 
P P 



sinh 



sinh 2i 

p 

xi - y' 



P 



(26) 



■^j;^P ■ In this way we get for 



tanh^ 

p 



l — xi 



tanh^ 
p J 



and altogether we get 

x-fF'ip) = 



n tanh ^ 
P p 



p2 ^ 



tanh ^ 
p p 



(27) 



(28) 



Note that the singularity present in the replicon part 
is canceled by the 1/p^ term in the longitudinal part, 
leaving a less divergent quantity in the long wavelength 
limit. Let us now look more in details at the infrared be- 
haviour of this correlation. In the far-infrared region, i.e., 
for momenta much smaller than the small mass 2xiwqi 
(see il2j), we have that xi/p ^ oo and 



p 



2xiwqi 



1/2 



1 

P 



(29) 



for <C 2xiwqi 



On the other hand, in the so called near-infrared region, 
i.e., in between the small mass 2xiwqi and the large mass 
2wqi, we have that xi/p <^ 1 and we get 



p2 -I- 2wqi{l — xi) 



(30) 



for 2xiwqi <C p^ <C 2xiwqi 



Thus, for intermediate momenta (or intermediate dis- 
tances, i.e., between the correlation lengths associated to 
the two masses), xJT'^ip) behaves as if it were massive. 
Only at very large distance does one feel the infrared 
singularity in 1/p (or the associated inverse power be- 
haviour in distance). The width of the region where the 
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apparent massive behaviour is felt depends on the size of 
X\. For 6 < d < 8, one has X\ ^ {w"^ /Cl t'^''-^^^/^, where t 
is the reduced temperature (Tc — T) /Tc and the num- 
ber of neighbours. As one reaches d — 6 from above xi 
remains smah hke w^/C- Below d = 6, xi is a function 
of w^/C with, at the critical fixed point, w^/C = e/2, 
i.e., xi = /(e) ~ e/2 for e <C 1. In three dimensions 
we have unfortunately no explicit knowledge of the size 
of xi, though the previous computation shows that, if 
the two mass-scales remain well separated, there exists a 
range of distances where massive-like behaviour occurs. 

In the far-infrared region one can apply scaling rules 
to guess what power-law behaviour arises below c? = 6. 
From the Kadanoff scaling hypothesis we expect the cor- 
relation function to decay as XfT^ip) = '4'{p O/P^'^i 
where ^ ~ is the associated correlation length. From 
Eq. ijSni) we then get 



Xi 



p 



that is 



p \2xiwqi 



1/2 




Xi p 



-V 



— (31) 



(32) 



or by Fourier Transform 



xr (0 



1 



(33) 



Consider now instead the near-infrared regime. In 
this region, the generalized correlation function with 
a = 2 shows an exponential decay with distance, the 
asymptotic power law decay arising only at still much 
larger length scales. As we shall see, this behaviour 
is peculiar of a = 2. Indeed for every other value of 
a the correlation function always decays algebraically 
(see Sec. IIII C 2(1 . Within the framework of the scal- 
ing/droplet model ^ .10] correlation functions always 
decay algebraically, with no difference between near and 
far-infrared. 

We note that in the near-infrared the a = 2 correla- 
tion function has a longitudinal nature in the usual field 
theory terminology. It is then appropriate to compare 
the present result with what happens in other systems 
where a continuous symmetry is broken. It is known that 
in 0{N) systems the longitudinal mode, which is mas- 
sive at the level of Gaussian fluctuations, is converted 
into a massless mode when loop corrections are taken 
into account Il7| | , and one may wonder whether the 
same also happens in spin glasses. The argument for the 
Heisenberg model, as detailed in Appendix relies on 
the fact that the effective interaction between Goldstonc 
modes vanishes in the infrared limit, and these modes are 
thus effectively free. In the Ising spin glass, however, the 
situation is different since in this case Goldstonc modes 
remain coupled in the infrared [T^ . As a consequence, 
the argument for the 0{N) case does not apply here and 



the longitudinal mode x"^^ in the near-infrared is ex- 
pected to remain massive: massive-like behaviour at in- 
termediate scales is expected to occur even if the loop 
corrections are included. 

Finally, we note that a different behavior is found if our 
analysis for d > 6 is extended to configurations which are 
constrained to have a fixed value of their overlaps. For 
example, for r = s = 0, that is for configurations having 
zero mutual overlap, the longitudinal part is identically 
zero, and we are left with a singular l/p^ behavior in the 
whole infrared region. (This happens in fact for all values 
of a, not just for a = 2.) 



2. The Case q / 2 

For a ^ 2 a calculation similar to the one detailed 
in the previous section enables us to estimate the most 
singular contributions in the limit p ^ 0. We find that 



a^2 



A -I- logp 



P- 



(34) 



where A and B are two numerical coefficients. These 
expressions are valid for d > 6. Using general scaling 
arguments for the integrals appearing in the calculation, 
at d < 6 one can argue that Eq. gives in real space 



a^2 



B 



(35) 



In d = 3, two re- 
= ±1 model give 



where i] is an anomalous dimension, 
cent numerical estimates for the Jtj 
ri = -0.35 ± 0.05 [3, and 77 = -0.22 ± 0.02 JJJ. The 
corresponding values for the exponent {i / 4){d — 2 + vj) of 
the leading term of Ea. (|35() are, respectively, 0.49 ± 0.04 
and 0.585 ±0.015. 



IV. COMPUTATIONAL METHODS 
A. Monte Carlo Simulations 

It is interesting to test the analytical predictions we 
have just provided (based on d > 6 calculations) to what 
happens in a three dimensional system. This means esti- 
mating the correlation functions of Eq. from numeri- 
cal simulations. To do that, we have generated 512 spin- 
glass samples, where the couplings Jy were randomly set 
to ±1; we considered mainly 12x12x12 cubic lattices 
with periodic boundary conditions (smaller lattices have 
also been used, allowing us to see the qualitative behavior 
of finite size effects). 

For each such sample we produced over 1000 equilib- 
rium spin configurations which were essentially indepen- 
dent (since they were separated by a large number of up- 
dates). This was achieved using the Parallel Tempering 
Monte Carlo method [23, |^ as follows. The tempera- 
tures used went from 0.1 to 2.0 in steps of AT = 0.1 (the 
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critical temperature Tc in this model is about 1.1). At 
each temperature we performed sweeps using the usual 
Metropolis algorithm and followed this by trial exchanges 
between neighboring temperatures; we call this sequence 
a pass of the algorithm. The value of AT was such that 
these exchanges succeeded with a probability close to 1/2. 

After thermalizing the system using a million such 
passes, we stored configurations every 1000 passes for fur- 
ther analysis. With these choices of parameters, we ver- 
ified that the configurations kept were compatible with 
being statistically independent by examining their over- 
laps. Furthermore, we also checked that each configu- 
ration visited many times all the temperature values as 
suggested in and that the distribution of overlaps for 
each sample was symmetric with high accuracy. Finally, 
the use of bit-packing enabled us to speed up the com- 
putation (this trick is particularly useful for the case of 
Jij = ±1 couplings). The overall computation was per- 
formed on a cluster of Linux-based PCs running at 333 
MHz and represents the equivalent of 2 years of mono- 
processor time. 



descent optimization method where we seek to have all 
mean overlaps be positive; the first pass is useful for gen- 
erating a good starting set. 

In many cases, we find that this iterative procedure 
converges to a fixed point where all mean overlaps are 
positive. This procedure has then split the configura- 
tion space into two parts that are related by the global 
spin flip symmetry. When the algorithm does not con- 
verge, we repeat the procedure with a new first pass after 
permuting randomly the configurations; this can lead to 
convergence, but if it does not, we simply stop after some 
number of trials and accept an "imperfect" decomposi- 
tion. 

Given the set of aligned configurations (about half of 
them have been flipped), we can then look at individual 
rather than mean overlaps. The main test of whether 
the configuration space decomposition is "perfect" is then 
the positivity of all the overlaps between these aligned 
configurations. 



B. Configuration Space Decomposition 



Parametrization of Correlation Functions 



All of the theory exposed in Section IIIII implies the 
presence of an infinitesimal magnetic field. Performing 
this limit numerically is a nuisance because of the neces- 
sity to have several values of the magnetic field and to 
control the accuracy of the limit. Fortunately, it is pos- 
sible to avoid such an extrapolation by remarking that 
the role of a very small magnetic field is to select half of 
the B = configurations. The simplest implementation 
of this is to keep only those configurations (generated at 
B = 0) whose total magnetization M is positive. This 
should give the correct selection when the lattice size L 
tends to oo, though of course with finite L this may lead 
to some undesirable effects. What properties does one 
expect in the thermodynamic limit? In the mean-field 
picture, there are many equilibrium states, all coming in 
pairs when i? = 0; if one adds an infinitesimal field, only 
one state in each pair survives and the overlap with the 
remaining states becomes positive. 

In practice, we have found that the selection of those 
configurations with M > left us with quite a few over- 
laps that were negative. (This was severe for some disor- 
der samples and less in others.) We thus sought a better 
decomposition of the configuration space into two parts 
related by the global spin flip symmetry. For each disor- 
der sample, we used the following procedure. (A nearly 
identical method was applied in 23]). 

In a first iteration, we successively consider the equi- 
librium configurations, and flip them if their mean over- 
lap with the previous configurations is negative. In the 
iterations thereafter, one repeats these flips but now us- 
ing the mean overlap with all the other configurations. 
In other words, we try to "align" the configurations as 
much as possible; effectively, our procedure is a kind of 



Within the theoretical analysis, one controls only the 
small momentum behavior of x" (c-g-? p ^ in for- 
mula (|34ll ). that is the large distance region. There is 
thus some arbitrariness in the choice of fitting functions 
when comparing the data to the theoretical predictions. 
Our choice of functions is as follows. Start with a cor- 
relation function G{p) on the infinite lattice that has a 
power law behavior p^"^ as p — > 0; the simplest choice is 



G(p) 



1 



2EJl-cos(p^) 



c/2 



(36) 



as motivated by the requirement that G{p) be periodic in 
the components of p. (We work in units of the lattice 
spacing.) This correlation function can be used also on a 
finite size lattice of size Lx Lx L with periodic boundary 
conditions as long as we take the correct discrete values 
of the Pf^ . We thus take this G{p) and Fourier transform 
it to coordinate space. All our measurements are for i 
and j on an axis of the lattice, and when performing the 
transform we must remove the p = contribution; this 
defines a function hereafter called Gc{r). An analogous 
procedure was used in j24j . One can also perform this 
computation with a mass term, leading to exponentially 
decaying correlation functions; wc do this by adding rri^ 
in the denominator of the right hand side of Eq. I|36p 
and setting c = 2, leading to the function Gmir)- (We 
performed checks on our code using the high precision 
values for such functions in psf.') The resulting family 
of functions we use for our fits are then a + bGc{r) and 
a + hGm{r). 
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FIG. 2: Example of distribution of overlaps for a disorder 
sample after the configuration space decomposition. T — 0.5. 

V. RESULTS 

A. Effect on Overlaps of the Decomposition 

In the thermodynamic hmit, we expect configuration 
space decomposition schemes to lead to aligned config- 
urations for which the overlaps are positive. However, 
one can expect some schemes to work better than others 
when L is relatively modest (in our case L = 12). As 
we already pointed out, the simplest procedure does not 
work so well, so let us examine the overlaps found with 
our "improved" scheme previously described. For each 
sample, we have viewed the distribution P{q) of overlaps 
found for each pair of configurations, both before and af- 
ter the alignment. In the majority of the samples, the 
procedure leads to nearly all overlaps being positive; yet 
there are some samples for which one still has a signifi- 
cant number of overlaps that are negative. In Fig.|2]we 
show an example where a rather small fraction of the 
overlaps are negative. When considering all 512 of our 
samples, the decomposition is reasonably good but not 
perfect; this is not surprising considering our lattice size 
(L = 12). Because of these effects, we expect our analysis 
of correlation functions to suffer from small systematic ef- 
fects since quantities such as (Si) are probably estimated 
with a residual bias. 



B. Correlation Functions 

Now we confront the different predictions of the replica 
field theory approach to the properties of the system as 
extracted from our simulations. 

The first important prediction of replica field theory is 
that the large distance limit of is (2 — a)/3 times the 
variance of q in the presence of an infinitesimal magnetic 
field. This prediction arises from the mean-field com- 
putation; it continues to hold when taking into account 
the Gaussian fluctuations; furthermore, as argued in the 
previous section, it should continue to hold within the 



0.03 - 
0.02 
0.01 - 



L , , , , , jJ 

0.4 0.8 1.2 1.6 2 2.4 

a 

FIG. 3: Large distance value 11 of the correlation function x" 
of Eq. ||7J as a function of a. Theoretically, the change of sign 
should be at a = 2. (The data are for T = 0.5.) 

loop expansion. We have fitted x"('i, j) using the class of 
functions a + hGc{r) and a -\- bGmir) defined before. We 
find that for a not close to 2 the power fit is far superior 
and gives a good value for x^- there is at least one or- 
der of magnitude difference between the power law best 
fit x'^ a-nd the exponential best fit x^- Oii the contrary, 
when a ~ 2, the exponential fit is better: in this case the 
decay is indeed so fast that only at distance 1 do we get 
sometimes a signal significantly different from zero over 
the plateau value. 

The constant value obtained from these fits gives the 
long distance limit of x"; for T = 0.5, we plot in Fig. O 
these limits as a function of a. The behavior is linear in 
a as it should, and the values change sign not far from 
a = 2. Furthermore, when plotting the data as a func- 
tion of (2 — a)/3, the slope should be (g^) — {q){q)- The 
numerical value we find for this slope is 0.040; this has 
to be compared to the numerically determined value of 
the variance of \q\ among our equilibrium configurations 
which is 0.047 at T = 0.5: these two values are numeri- 
cally close. Furthermore, since we find that this variance 
decreases as L increases, the two quantities probably do 
become equal at large L. Our measurements thus sup- 
port the theoretical claim that the large distance limit of 
X" is exactly given by Eq. H21|l and are consistent with 
a large body of published results "^l. The situation is 
similar at other temperatures in the low T phase. Notice 
that the plateau value at a = 2 is different from zero: 
the effect is small, and on the L — 12 data we get a zero 
plateau close to a = 2.4. But this discrepancy decreases 
with increasing lattice sizes, making manifest its finite 
size nature. 

Let us move on to the next prediction of the analytic 
computation concerning how the correlation functions 
tend toward their large r limit. In Fig. ^ we show the 
raw data for x^C*^) for a = 0, . . . , 2 in steps of 0.5. We 
have also displayed the fits obtained when using the func- 
tional form a + bGc{r) for the data 1 < r < 6. For these 
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FIG. 4: Data for x"(r) for a = 0, . . . , 2 in steps of 0.5 (from 
top to bottom). For a ^ 2 we display the best fits using 
a + bGc{r). (The data are for T = 0.5.) 

values of a, the best values of the power c are 2.51, 2.49, 
2.45 and 2.34 for the first four. These fits are good, in line 
with the absence of visible systematic errors. On the con- 
trary, when using the exponential fits, i.e., a + hGm(,r), 
the fits are not very good when a is away from 2; in 
fact, the is minimized when the mass tends to zero 
for a < 1.5. Finally, when a — 2, the exponential fit is 
good and leads to a mass of 1.1. 

Put all together, these analyses strongly favor a power 
law decrease with r when a ^ 2 with a power that is 
consistent with being a-independent. To leading orders 
one then has 

x"^^«-^-f,| + ^ (37) 

where crj^^| is the variance of \q\. In our data, A = c? — c is 
compatible with being an a-independent constant whose 
value is close to 0.5. This value may be compared to the 
value 9 = 0.20 expected within the scaling/droplet pic- 
ture, and to the exponent (3/4) (d— 2 + 77) predicted using 
scaling arguments within the replica field theory calcula- 
tion (equal to 0.49 or 0.585, according to the numerical 
estimate of 77 used in the formula - see Sec. IIII C 2|l . In 
contrast, when a = 2, our data favor instead the fit to 
an exponential decrease with r, and thus supports the 
extrapolation of the d > 6 theoretical analysis of the 
near-infrared behaviour down to d = 3. 

VI. DISCUSSION 

We have investigated the properties of certain correla- 
tion functions in spin glasses involving two spins. As a 
first result, we have found that the large distance behav- 
ior of the two connected correlation functions Gi and G2 
(see Eq. Q) satisfies 

aGi(oo) + (l-a)G2(oo) = ^ [F> - ¥'] (38) 



and thus Gi(cx)) = o'j^^|/3 and 62(00) = 2CTj^^|/3 where 
aj^^l is the variance of the absolute value of the overlap, 
\q\. (In all our computations, we assumed the presence 
of an infinitesimal field: the net effect of that field is to 
decompose the configuration space, rendering overlaps 
positive.) 

A second series of results concerns the way the corre- 
lation functions Gi and G2 tend toward their limits. We 
found numerically a decrease close to 1/r^/^ which is a 
bit different from the droplet/scaling prediction of 1/r^ 
since 9 « 0.2. This power law decay disappears when 
considering precisely a = 2, i.e., 2Gi(r) — G2(r): in that 
case one has both a limiting value compatible with zero 
and a very fast decay, suggestive of an exponential law. 

All of our analytical predictions are based on replica 
field theory calculations using an effective Lagrangian; 
this approach should be reliable in dimensions d > Q. 
Nevertheless, the main predictions, namely the limiting 
values of Gi and G2 and their approach to their limits, 
all are corroborated in d = 3 by our numerical study. 
Also the prediction of a massive-like regime for 2Gi — G2 
seems to be nicely consistent with the numerical findings 
that show for this correlation a very fast decay. 

It is of interest to note that Fisher and Huse 0, 0| 
and Bray and Moore had also suggested that there 
could be a peculiar behaviour for a longitudinal-like lin- 
ear combination of Gi and G2, due to the cancellation 
of the leading term. More precisely, the claim there was 
that the power law behaviour would be changed from 
6* « 0.2 to a larger value (determined by the expansion 
of the zero temperature distribution of the internal fields 
around h = Q - see Sec. 4.1 of 0)- We do not find a 
decrease compatible with 1 /r^ for Gi and G2 separately, 
while we find a decay compatible with an exponential law 
for 2Gi — G2, at least in the near infrared; thus qualita- 
tively, the correspondence with the droplet is not good. 

Coming back to our results for the large distance lim- 
iting values of Gi and G2, we can derive a simple "sum 
rule" as follows. Each term in x"(r) can be computed 
when |j — j| — > 00 via replica moments just as we did in 
Eq. P: 



^Ihn {S,S,){SCi{S,) = + (q)2] ^ (q.^^) , (39) 

and 

hm {S,)HS,Y = ^[2Fy + ¥^] - , (40) 

where, in the last equalities, the averages (• • •) are per- 
formed with respect to a replicated equilibrium measure. 
We then obtain the following sum rule when using Eg. 1381 
at a = 2: 



(g22)-4((?i2gi3>+3(gi2-734>=0. (41) 

It is also possible to derive this relation using arguments 
based on stochastic stability 
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As a final comment, note that within the scahng or 
droplet pictures, when the distance between i and j di- 
verges one has 



lim 



{s.y {s,y^ 



(42) 



This does not hold in the presence of replica symmetry 
breaking, and instead we have a relation following from 
Eq. pn|). 

Several questions remain open. One would like to de- 
termine analytically the power law decrease of Gi and 
G2 which here was compatible with 1/r^/^; is that the 
exact value, and what is it in higher dimensions? An- 
other issue is related to the presence of the massive-like 
behaviour for the a = 2 correlation function in the near- 
infrared. It would be important to get some analytical 
estimates of the range where this behaviour does hold 
also in d = 3. Finally, one may wonder whether there is 
any deep reason for the cancellation of the leading parts 
of Gi and G2 for a = 2. One may imagine that it is 
somehow related to sum rules. Then one may conjecture 
that to each sum rule for overlaps (obtained for instance 
using stochastic stability), there is an associated corre- 
lation function which decays to zero anomalously fast. 
This conjecture should be testable for several sum rules 
using Monte Carlo methods. 
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APPENDIX A: REPLICA FOURIER 
TRANSFORM 



The Replica Fourier transform is a powerful technique 
which greatly simplifies the diagonalization and inver- 
sion of the fluctuation matrix G (see Eq. (|10|) ')'). The 
idea is the following. Convolutions in configuration space 
become products in Fourier space when appropriately 
defining a Fourier transform in replica space; thus sums 
over replica co-distances may be transformed into associ- 
ated products. Similarly, matrices can be inverted in the 
replica Fourier space and then transformed back to get 
the desired expression. 



Let us assume that replica objects (matrices, tensors 
etc.) only depend on replica co-distances, and belong 
to the hierarchical structure described in Sec. IIIII The 
Replica Fourier Transform (RFT) is a discretized 
and generalized version of the algebra introduced in [23| . 
For objects like qab which depend on a single replica co- 
distance aOb — t the RFT is defined as 



_R+1 
t=fe 



qt-i) 



(Al) 



where the pt are the sizes of the Parisi blocks, and objects 
with indices out of the range of definition are taken as 
equal to zero. 

The inverse transform is given by: 



Qt 



t 

E 

fc=0 



1 

Pk 



ilk 



(A2) 



With these definitions, the convolution J2c ^acBcb = Cab 



becomes in RFT Af,Bf, 



C^. However, in our case the 



problem is slightly more complicated. Starting from the 
known expression of the fluctuation matrix M of Eq. I|18|l . 
we want to compute the propagator G = M^^ . Since 
both M and G are tensors bearing four indices, the uni- 
tarity equation f^fab;cdQcd;ef _ involves a dou- 

ble convolution in replica space. If the fluctuation matrix 
and the propagator are expressed in terms of replica co- 
distances, as described in Eqs. ((T^ and (^0)) . one needs 
to Replica Fourier transform with respect to lower indices 
(i.e., the indices referring to the cross-overlaps between 
pairs of replicas). 

In the Replicon sector, the double transform reads: 

R+l R+l 

Mu = (A3) 

u—k v—l 

with, by definition, fc, / > r + 1. The inverse transforma- 
tion can be obtained by applying twice the inverse RFT 
defined in Eq. ljA2|l . From the Lagrangian, Eq. ©, we 
can easily recover the explicit expression of the fiuctua- 
tion matrix. One has (with the notation of Sec. IIIIB|l : 



M 



R+1.B.+1 



p'-2iT + uqi) 
-wqu , 
~wqy . 



(A4) 



From this, by noting that the mean-field solution gives 
in Replica Fourier space t = — w^/g, one gets: 



^'^U + ^[(96 - 9fe) + (96 - 9/)] - 2ug^ 



(A5) 



In the limit of an infinite number of steps of replica 
symmetry breaking, where now r e [0,1] 

while q{r) obeys Eq. (|17|l . Correspondingly, the Replica 
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Fourier transform becomes g^, 
and 



J^^ dxxdq/dx — 2uq{xi) 



G^P) 



71 r' ^ 



4u 



[k^ + 1^ - 2r^ 



(A6) 



Finally, by simple inversion, we get the replicon propa- 
gator of Eq. 

The computation of the longitudinal Fourier com- 
ponent of the propagator is less straightforward. In 
the longitudinal-anomalous sector the fluctuation matrix 
Mp"^ carries only a lower index (cross-overlap), and thus 
only one Replica Fourier Transform is needed. However, 
the lower index runs over a hierarchical tree and when t 
crosses the upper 'passive' indices, the multiplicity may 
change (sums over replica co-distances such as r, s,t are 
equivalent to the original sums over replica indices only 
if the correct multiplicity is taken into account). To deal 
with this, one has to generalize the definition of Replica 
Fourier transform in the presence of passive indices: 



k 



(A7) 



t=k 



with, for r < s. 



(r,s) 
Pt ^ Pt 



t<r, 
2pt r <t < s , 
Apt r < s < t . 



(A8) 



With such a definition, the unitarity equation in the 
longitudinal-anomalous sector becomes: 



k 



Z — / k 4 

t=0 



k > 



(A9) 



where 



and 



gl = G!^I k<r + l, 

^k r+l,r+l — ' 

= G!^ ^ fc > r + 1 , 

r+l,k ' 



5'^ 



t -pf-pf+l 



(AlO) 



(All) 



This equation can be solved with the use of Gegcnbauer 
functions whenever M^'^ depends only on min(r, s) or on 
max(r, s) (see ^3 for details). The result for G^^^ is the 
one reported in Eq. H2()(l . 

Finally we note that, in the replicon geometry (r = s), 
there is a Longitudinal- Anomalous (LA) contribution to 
the replicon fiuctuation matrix and to the propagator. If 
we write 



~-R +A , 

the LA part can be shown to be 14] 



(A12) 



m: 



(A13) 

When performing a double RFT as in Eq. (EU, the LA 
contribution disappears since it bears only one lower in- 
dex and Eq. I|A3|I is thus a purely replicon contribution. 




FIG. 5: Series expansion for the longitudinal propagator in 
the Heisenberg model. 



APPENDIX B: THE CASE OF THE 
HEISENBERG MODEL 

Let us consider the case of the isotropic Heisenberg 
model with N components. The Hamiltonian for this 
system reads: 



<5 



[ i=i ■ \i=i 

(Bl) 

In the low temperature region, for d > 2, this system 
develops a spontaneous magnetization M . By expanding 
Eq. ljBl|) around the mean-field solution AP — 6\fi\/g, 
i.e. ^ = M + <j), we get a field theory where fiuctuations 
along the direction of M (i.e. longitudinal modes) are 
massive, while fluctuations transverse to it are massless 
(zero or Goldstone modes). More precisely, 



n = n 



mf 



(V0t)' 



(B2) 



where 0l is the longitudinal mode, and 0t the (N — 1)- 
dimensional transverse mode. The free propagators then 
read: 



Glip) 

gUp) 



fM2 



(B3) 



We now show with a simple argument that, when con- 
sidering loop corrections, the longitudinal mode, which 
is massive at the bare level, becomes massless. To see 
that, let us exhibit the series expansion of the longitudi- 
nal propagator Gl as in Fig. [S] Here the L- lines corre- 
spond to G'j^ and the wavy lines stand for the magneti- 
zation M . Since we are interested in the p — > limit of 
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r(p) - 

T 




FIG. 6: Effective interaction between transverse modes, ob- 
tained by integrating out tlie longitudinal modes. 

Gl(p), all other lines are transverse lines, G^. The thin 
vertices are couplings g, the heavy vertices correspond to 
the effective coupling Ft between transverse modes. We 
get, integrating out the transverse loops, 

(B4) 

where / is a numerical prefactor. The effective transverse 



interaction can be obtained by integrating out the con- 
tribution of the longitudinal modes (see Fig. ^ , yielding 



^ -9 - YWTW ^ P^^W ■ ^^^^ 

Thus, TTij}) vanishes in the infrared limit p — > and 
the Goldstone modes are effectively free. Inserting this 
expression in Eq. l|B4p . we see that the massive bare be- 
havior l/[p^ -I- (5/3)M^] of the longitudinal propagator 
is changed into a massless one l/p'^"'^. 

This simple argument cannot be applied to the Ising 
spin glass. Indeed there the Goldstone modes are the 
bottom of a (transverse) band with no gap separating 
the massless from the massive modes (see As a 

result, zero modes remain coupled in the infrared, i.e., 
^t{p = 0) 7^ 0, and they develop an anomaly. 
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